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Abstract 


The  class  of  time-varying  linear  systems  which  are  two-time-scale  on 
an  interval  may  be  decoupled  by  a  time-varying  transformation  of 
variables  into  separate  subsystems  containing  the  slow  and  fast  dynamic 
parts.  The  transformation  is  obtained  by  solving  a  nonsymmetric  Riccati 
differential  equation  forward  in  time  and  a  linear  matrix  differential 
equation  backward  in  time.  Small  parameters  are  identified  which  measure 
the  strength  of  the  time  scale  separation  and  the  stability  of  the  fast 
subsystem.  As  these  parameters  go  to  zero,  the  order  of  the  system  is 
reduced  and  a  useful  approximate  solution  to  the  original  system  is 
obtained.  The  transformation  is  illustrated  for  examples  with  strong  and 
weak  fast  subsystem  stability. 
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1 .  Motivating  examples  of  time-invariant  problems. 

The  longitudinal  dynamics  of  an  F-8  aircraft  (cf.  Etkin  [12]  and 

Teneketzis  and  Sandell  [35])  can  be  modeled  by  an  initial  value  problem 

for  a  fourth  order  linear  system  of  the  form  x  =  Ax  +  Bu  with  the 

physical  variables  being  the  "primarily  slow"  velocity  variation  and 

flight  path  angle  and  the  "primarily  fast"  angle  of  attack  and  pitch 

rate.  The  single  control  is  the  elevator  deflection.  The  exact  solution 

for  the  free  response  of  the  components  of  x(t)  is  plotted  in  Figures  1 

and  2.  Our  objective  is  to  determine  a  solution  x(t)  of  a  reduced  second 

order  model  which  will  approximate  the  dynamics  of  the  given  fourth  order 

model  away  from  the  initial  time  t  =  0.  We  wish,  in  particular,  to  avoid 

integration  of  the  full  order  system  or  a  complete  eigenanalysis  of  A 

At  A(  t-s ) 

which  would  provide  the  exact  solution  x(t)  =  e  x(0)  +  e  v  'B(s)u(s)ds. 

J0 

At 

We  note  that  approximations  to  the  matrix  exponential  e  are,  indeed, 
still  the  subject  of  current  research  (cf.  Moler  and  Van  Loan  [23])  and 
that  they  are  not  simple  to  compute. 

Our  criteria  for  such  approximations  will  naturally  involve  the 


eigenvalues  of  A.  For  this  problem,  we  have  the  "slow"  eigenvalues 
s^  2  =  -0.0075  +  i 0 . 076  and  the  "fast"  eigenvalues  f]  ^  =  -0.94  +  i 3 . 0 . 
Our  method  will  rely  upon  the  time-scale  separation,  measured  by  the 
smallness  of  the  parameter  p  =  =  0.024,  and  the  relative 

stability  of  the  fast  eigenvalues,  measured  by  the  parameter 
=  - [ Re  S2 | /Re  f-j  =  0.0081.  Most  important,  however,  is  that  the  ratio 
of  the  fast-mode  decay  constant  ( - 1 / Re ( f ^ ) )  to  the  length  T  of  the  time 
interval  of  interest  satisfies  cr  .?  T "To  tw*TV  ^  « 
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For  p  and  e  small  we  expect  our  reduced  order  model  to  be  a  good  approxi¬ 
mation  to  the  solution  on  an  interval  0(eT)  <  t  <  T,  while  on 
0  <  t  <  0(eT)  any  fast  mode  components  excited  by  the  initial  conditions 
may  be  significant  and  the  approximation  x(t)  which  ignores  them  would  be 
inappropriate.  As  the  figures  suggest,  T  must  be  quite  large  in  order 
for  the  initial  layer  to  be  relatively  narrow. 

For  large  dimensional  linear  problems,  one  cannot  readily  compute 
exact  solutions  to  compare  approximate  solutions  against.  In  power 
system  models,  systems  involving  several  hundred  variables  are  common. 
They  are  often  approximated  by  reduced  order  models  involving  both 
differential  and  nonlinear  algebraic  equations  which  neglect  fast  initial 
transients  (cf.  Van  Ness  [37]).  An  algebraic  system  g(x,z,t)  =  0  could 
correspond  to  a  steady-state  for  the  differential  system  cz  =  g(x,z,t) 
with  the  small  parameter  e  representing  "parasitics".  The  practical 
importance  of  obtaining  reduced  order  models  follows  largely  because  the 
computational  effort  involved  in  numerically  integrating  systems  of 
differential  equations  increases  at  least  as  the  square  of  the  order. 

A  second  example  of  a  two-time-scale  problem  is  the  sixteenth  order 
model  of  a  turbofan  engine  which  was  the  theme  problem  for  a  recent 
conference  on  control  of  linear  multi-variable  systems  (cf.  Sain  [31], 
Skira  and  DeHoff  [33],  and  DeHoff  and  Hall  [10]).  The  linear  model  is  of 
the  form  x  =  Ax  +  Bu,  y  =  Cx  +  Du,  with  the  state  variables  being  fan 
speeds,  pressures  and  temperatures.  The  five  control  variables  u  are 


%\ 


4 


fuel  flaw,  exit  nozzle  area,  two  vane  angles,  and  compressor  bleed;  and 
the  outputs  y  are  thrust,  total  airflow,  a  temperature,  and  two  stall 
margins.  The  objective  in  the  controller  design  is  to  achieve  rapid 
thrust  response  without  violating  several  operating  constraints.  The 
model  is  one  of  thirty-six  different  linear  models  obtained  from  a  non¬ 
linear  simulation  of  the  engine.  It  represents  the  turbofan  operating  at 
sea  level  with  near  maximum  non-afterburner  power.  This  is  an  operating 
point  which  every  engine  must  pass  through  at  takeoff.  Based  on  the 
eigenvalues  of  this  model  and  T  =  2  the  time-scale  separation  and  fast 
mode  stability  parameters  are  (u,e)  =  (0.304,  0.000867),  (0.371,  0.0285), 
(0.383,  0.0744),  respecti vely,  for  the  number  n^  of  "slow"  modes  chosen 
as  15,  5  or  3.  Since  an  order  reduction  from  16  to  15  isn't  substantial, 
we  shall  use  =  5.  In  all  cases,  the  time  scale  separation  and  relative 

stability  parameters  p  and  o  are  only  marginally  small,  while  the  fast-mode 
stability  parameter  e  is  quite  small.  We  nonetheless  obtain  good 
approximate  solutions  by  solving  a  reduced  (fifth)  order  system  instead 
of  the  original  sixteenth  order  problem.  In  Figures  3  and  4,  the  exact 
solution  of  the  sixteenth  order  problem  and  the  solution  of  our  fifth 
order  model  of  the  slow  dynamics  are  plotted  for  the  thrust  and  fan  speed 
in  response  to  changes  in  controls.  The  control  inputs  are  step  changes 
in  fuel  flow  rates  and  inlet  guide  vane  position.  The  second  case,  cf. 
Figure  4,  provides  a  severe  test  to  the  reduced  order  model  since  the 
inlet  guide  vane  is  located  at  the  front  end  of  the  engine  and  there  is 
some  delay  before  its  effect  is  propagated  to  the  net  thrust.  We  note 
that  the  approximations  are  not  good  for  t  <  0.28  =  10c  and  that  this 
initial  layer  will  become  narrower  relative  to  T  as  c  tends  toward  zero. 
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Figure  3: 

Turbofan  engine:  Respoi 
of  thrust  and  fan  speed 
a  500  Ib/hr  step  change 

fuel  flow  rate.  | 
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Figure  4: 

Turbofan  engine:  Respor 
of  thrust  and  fan  speed 
a  10  degree  step  change 
in  inlet  guide  vane 
position. 
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2 .  The  time-varying  problem— an  exact  approach. 

Several  earlier  papers  (cf.  Kokotovic  [16],  Chow  and  Kokotovic  [7], 
O'Malley  and  Anderson  [28]  and  Anderson  [3])  have  discussed  time- 
invariant  problems,  so  let  us  now  consider  the  time-varying  system 

(1)  x  =  A(t)x  +  B(t)u(t)  ,  0  <  t  <  T  , 

where  A  and  Bu  are  specified. 

A  system  such  as  (1)  will  be  called  two-time-scale  on  the  interval 
[0,T]  if  the  spectrum  X(A(t))  of  the  n  x  n  matrix  A  can  be  partitioned 
into  two  sets  S(t)  and  F(t)  with  and  n^  =  n  -  n.j  elements, 
respectively,  such  that 

X ( A ( t ) )  =  S(t)  U  F(t )  , 
where  the  eigenvalues  satisfy 

(2)  max  |s,(t)i  =  s(t)  «  f(t)  =  min  jf.(t)j 

Siss  ffF 

throughout  0  <  t  <  T  with 

(3)  y  =  max  (s(t)/f(t))  <<  1  . 

0<t<T 

Roughly,  then,  p  is  an  upper  bound  for  a  ratio  of  time-varying  eigen¬ 
values.  We  note  that  if  j Re  f.(t)|  is  large,  a  corresponding  vector 

0 

solution  of  the  unforced  system  will  be  locally  exponential ly  growing  or 
decaying,  while  if  j Im  f j ( t ) |  is  large,  the  corresponding  solution  will 
oscillate  rapidly  locally.  We  also  note  that  different  modelers  might 
select  different  values  of  n-j  for  the  same  system,  and  that  the  more 
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difficult  problems  where  varies  across  [0,T]  will  not  be  considered 
here.  Finally,  the  common  situation  where  the  eigenvalues  of  A  cluster 
into  several  sets  might  be  handled  by  repeated  application  of  our 
technique,  cf.  Kokotovic  et.  al_.  [17]  and  Winkelman  et  al.  [40]. 

For  general  time-varying  systems,  it  is  well  known  that  eigenvalue 
stability  does  not  imply  stability,  cf.  e.g.  Coppel  [8],  The  result  is, 
however,  more  nearly  true  for  singularly  perturbed  systems.  Thus,  for 
the  singularly  perturbed  initial  value  problem  for 

x  =  A(t,(c)x  +  B(t,K)z  +  C(t,i<)  , 

kz  =  D(t,«)x  +  E(t,ic)z  +  F(t,<)  , 

with  smooth  coefficients  on  0  £  t  £  T,  the  limiting  solution  as  <  tends 

to  0+  on  an  interval  0  <  t  £  T  will  satisfy  the  reduced  order  system 

X  =  A(t,0)X  +  B(t,0)Z  +  C(t ,0)  ,  X(0)  =  x ( 0 )  , 

0  =  D(t ,0)X  +  E(t,0)Z  +  F(t,0)  , 

provided  the  matrix  E(t,0)  is  stable  throughout  0  £  t  <_  T.  Further,  an 
initial  boundary  layer  (or  region  of  nonuniform  convergence)  occurs  in 
the  z  variable  near  t  =  0  and  the  fast  dynamics  there  evolve  on  a 
t  =  t/<  time  scale,  cf.  O'Malley  [26,27].  Such  theory  suggests  that 
eigenvalue  stability  may  be  appropriate  for  determining  the  behavior  of 
two-time-scale  systems.  These  results  apply  to  systems  where  the 
coefficient  matrices  A,...,F  have  bounded  t  and  k  derivatives.  Related 
problems  on  the  semi -inf ini te  interval  t  >  0  are  treated  in  Hoppensteadt 
[15],  Barman  [5],  and  Vidyasagar  [36].  With  less  smoothness. 
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counterexamples  exist  and  caution  must  be  observed,  cf.  Kreiss  [18,19].  For 
these  reasons,  Kreiss  introduced  hypotheses  demanding  that  E(t,0)  be 
"essentially  diagonally  dominant." 

We  shall  not  suppose  that  the  given  system  (1)  is  two-time-scale, 
but  rather  that  it  can  be  transformed  into  a  two-time-scale  system  by  a 
time-varying  transformation 

(4)  y  =  T(t)x  , 


with  the  system  for  y  being  "time-scale  decoupled"  throughout  0  £  t  £  T. 
Specifically,  let  the  transformation  matrix  T  have  the  form 


(5) 


T(t) 


In  +  K(t)L(t) 
nl 

L(t) 


K(tT] 

I 

! 

1  0  1 
"l 

I 

= 

0  1 

■ 

L  ( t )  I 

n2_ 

L  n2j 

1 

,  1 

n 

2J 

and  let  the  matrices  L(t)  and  K(t)  be  determined  so  that 


(6)  y  =  A(t)y  +  T(t)B(t)u  , 


where  A  has  the  block-diagonal  form 


A(t) 


An(t)  0 
0  A22(t) 


with  the  n-j  eigenvalues  of  A^  being  small  in  magnitude  compared  to  those 
of  A^2  throughout  0  <  t  <  T.  For  u  =  0,  the  slow  modes  for  (6)  would  be 


given  by  y  = 


where  y^  satisfies  the  lower  order  system  y^  =  A^y^ 


J 
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while  the  fast  modes  are 


ly2J 


where  =  A^y,,.  We  note  that  the  trans¬ 


formation  matrix  T  has  the  explicit  inverse 


(7) 


-L(t) 


-K(t) 

I  +  L(t)K(t) 
n2 


so  T  is  always  nonsingular  and  transformations  between  x  and  y  coordinates 
are  particularly  convenient.  Analogous  transformations  have  been  employed 
in  the  singular  perturbati ons  context  by  Wasow  [38],  Harris  [14],  and 
Kokotovic  [16],  for  discrete  problems  by  Phillips  [29],  and  for  difference 
equations  by  Matheij  [21]. 

As  a  first  step  toward  time-scale  decoupling,  let  us  set 


(8)  z  =  T-j  (t)x 


for  the  block  triangular  matrix 


^(t) 


L<‘> 


Clearly 

(9)  z  =  A ( t ) z  +  T-|  (t  )B(t  )u 


where 


9 

A(t)  =  (T -| A  +  T-j  )T^ 1  -  (A.j) 

An  "  A12L  A12 

L  +  LA-j  ^  -  A^L  -  LA^L  +  A^-j  A 22  +  LA^ 

presuming  the  original  A  matrix  and  A  are  both  partitioned  after  their 
first  n,  rows  and  columns.  In  order  for  A  to  be  upper  block-triangular, 
the  r*2  x  n^  matrix  L  must  satisfy  the  matrix  Riccati  equation 

(10)  L  =  A22L  -  LAn  +  LA]2L  -  A21 

throughout  0  £  t  £  T.  Selecting  L(e)  =  0  for  a  yet-unspecified  endpoint 
e  =  0  or  T  makes  ^(e)  =  0  and  T  (e)  a  similarity  transformation.  Thus 
A(e)  will  be  two-time  scale  provided  A(e)  is.  Let  us  suppose 

(HI)  A(e)  has  n^  “slow"  eigenvalues  in  S(e)  and  n2  "fast" 
eigenvalues  in  F(e),  e  =  0  or  T. 

This  will  actually  determine  the  integers  n]  and  n2  used  throughout.  We'll 
later  find  that  selecting  e  r  0  (e  =  T)  will  be  natural  if  the  fast  eigen¬ 
values  of  A(t)  are  all  stable  (unstable)  everywhere. 

We  now  begin  an  extended  discussion  on  how  to  compute  L(e).  In  so 
doing,  we  make  improvements  on  previous  solutions  to  the  time-invariant  A 
problem  for  which  L(t)  is  constant.  If  we  partition  the  spectral  decompo- 
si tion  of  A(e)  as 

P,  °1  , 

A(e>  =  Mi_0  ' 

where  '.(J-j)  =  S(e)  and  M  =  (M^),  we  can  always  reorder  the  entries  in 
the  state  vector  x  so  that  the  n^  *  n^  matrix  is  nonsingular.  The 


I  u 


pll 

columns  of  L  will  span  the  n,  dimensional  eigenspace  of  A(e) 

|_  21 J  1 

corresponding  to  the  slow  eigenvalues  in  S(e)  and 

(11)  L(e)  = 

will  be  the  unique  solution  of  the  algebraic  Riccati  equation 

(12)  A22(e)L  -  LAn(e)  +  LA12(e)L  -  A2](e)  =  0  , 
i.e.  L(e)  =  0,  achieving  the  time-scale  decoupling 

(13)  A(An(e))  =  S(e)  and  \(A22(e))  =  F(e)  . 

Though  the  matrix  equation  (12)  has  many  solutions,  only  (11)  provides 
the  desired  time-scale  separation  (cf.  Anderson  [2]).  We  also  note  that 
(11)  avoids  the  use  of  vectors  in  the  n2  dimensional  fast  eigenspace. 

An  alternative  representation 

(14)  L(e)  =  Q‘^Q21 

in  terms  of  the  left  eigenspace  corresponding  to  the  n2  fast  eigenvalues 
of  A(e)  would  be  more  practical  if  n2  «  n-j .  The  corresponding 

upper  triangular  transformation  might  then  be  more  convenient  then  , 

since  it  would  first  isolate  the  purely  slow  component.  Here  we  have  parti¬ 
tioned  M  ^  =  (Q^j)  after  its  first  n-j  rows  and  columns  and  the  inverti  bi  li ty 
of  Q22  follows  from  that  of  H^.  The  nontrivial  result  (11)  follows  via 
linear  algebra,  as  does  (14).  Specifically,  if  A^(e)  has  the  decompos¬ 
ition  XGX  \  the  algebraic  Riccati  equation  can  be  rewritten  as 


„  -/ 


tv 


X 

Y 


A21(e)  -  A99(e)L(e)  =  -L(e)XGX_1.  For  Y  =  -L(e)X,  A(e) 


22 


A(A-j-j(e))  =S(e)  implies  that  we  must  have 


"11 

LM21. 


G,  so 


K  for  some  non¬ 


singular  K,  i.e.  L(e)  =  -YX  ^  Calculating  further  with  this 

L(e),  A(e)  =  T-|(e)A(e)T^  (e)  is  upper  block  triangular.  Recent  work  in 
Medanic  [22]  also  describes  the  invariant  manifolds  of  such  matrix 
Riccati  equations.  Watkins  [39]  mentions  numerical  difficulties  occurring 
when  M,-j  is  ill-conditioned. 

Note  further  that  any  n-j  dimensional  basis  of  the  slow  eigenspace 

Pm, 


could  be  used  instead  of 


_M21-i 


in  (11)  to  obtain  L(e).  One  possibility 


is  to  compute  n^  Schur  vectors  for  this  slow  eigenspace,  cf.  e.g.  Laub 
[20].  Once  an  approximate  L(e)  is  obtained,  one  may  improve  on  its 
accuracy  by  iteration  in  the  linear  equation 


(15) 


-1, 


Li+1  =  (A22(e)  +  LiAT2(e) )~  (LiAl 1 (e)  +  A2l(e)) 


Anderson  [2]  shows  that  this  iteration  converges  linearly  with  asymptotic 
rate  f  \e)~s(e) ,  so  this  method  is  particularly  well-suited  to  systems 
whose  time-scale  separation  parameter  u  is  very  small. 

The  iteration  scheme  (15)  can  be  obtained  from  the  simultaneous 
iteration  method  (cf.  Stewart  [34]  and  Avramovic  [4])  for  calculating 
the  dominant  eigenspace  corresponding  to  the  n^  fast  eigenvalues  of  A(e). 
That  method  generates  the  n2  x  n  matrix  V  as  the  limit  of  the  iteration 


(16)  Vk+,  =  VkA(e)  . 
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Splitting  Vk  =  [V^  Vk2]  after  its  first  n-j  columns  and  setting  = 
V^Vkl  (cf.  the  alternative  representation  (14)  for  L(e)),  (16)  reduces 
to  (15).  The  asymptotic  rate  of  convergence  f~^(e)s(e)  was  known  in  this 
context.  Indeed,  the  fact  that  (16)  converges  globally  under  very  mild 
assumptions  on  VQ  implies  that  our  iteration  scheme  (15)  will  also  be 
robust  with  respect  to  initial  iterates.  Thus  Lq  need  not  be  generated 
through  a  preliminary  approximate  eigenanalysi s  of  the  slow  eigenspace 
of  A(e);  in  practice  a  trivial  LQ  achieves  convergence.  The  reader  won't 
be  surprised  to  find  closely  related  analysis  in  the  stiff  differential 
equations  literature,  cf.  e.g.  Alfeld  and  Lambert  [1]. 

The  Riccati  differential  equation  for  L(t)  will  have  the  constant 
solution  L(e)  when  A  is  constant  or  when  it  is  possible  to  find  a 


constant  basis 


with  N-j  invertible,  for  the  slow  eigenspace  of  A(t). 


Otherwise,  we  need  to  integrate  the  n^  *  n^  dimensional  initial  or 
terminal  value  problem  (10),  (11)  for  L(t).  We  shall  assume  that  it 
provides  a  transformed  system  for  z  which  is  two-time-scale.  Specifically, 
we  suppose: 


(H2)  the  solution  L(t)  of  the  matrix  Riccati  problem  remains 
bounded  throughout  0  <_  t  £  T  and  the  eigenvalues  of  the 
matrix  A^(t)  =  A^  -  A^L  remain  small  in  magnitude 
compared  to  those  of  A00(t)  A22  +  LA.J2  throughout  the 

interval . 


» 
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If  this  hypothesis  fails  at  any  point,  our  order  reduction  procedure 
will  simply  not  be  appropriate.  We  note  that  some  leeway  is  allowed  in 
judging  the  separation  of  eigenvalues  between  and  »  i.e* 
deciding  just  how  small  a  p  is  small  enough.  Computational  and  stability 
aspects  of  the  integration  procedure  will  be  illustrated  below  through 
discussion  and  examples. 

One  can  proceed  further  and  block  diagonalize  the  upper  triangular 

A 

matrix  A  by  a  second  nonsingular  transformation 

(17)  y  =  T2{  t)z  =  T(  t)x 

for 

T( t)  -  T2(t)T}(t) 
and 


cf.  (5).  Thus  (17)  converts  (1)  into  the  two-time-scale  system 
(18)  y  *  A(t)y  +  B(t)u 


cf.  (6)  where 


A(t)  =  (TA  +  T)T_1 


An(t) 

0 


Alz(tf 

A22(t) 
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with  =  K  -  A-j -j K  +  KA22  +  and  B  =  TB.  If  the  n-j  *  matrix  K 
satisfies  the  linear  differential  equation 

(19)  K  =  An(t)K  -  KA22(t)  -  A12(t)  , 

the  matrix  A  will  be  block  diagonalized  and  the  system  for  y  will  be  time- 
scale  decoupled,  i.e.  the  system  for  the  first  n-j  components  of  y  will  be 
completely  decoupled  from  that  for  its  last  n2  components.  Corresponding 
to  the  endpoint  condition  L(e)  =  0  for  L,  we  now  impose  the  condition 
K(T  -  e)  =  0  at  the  opposite  endpoint  because  the  variational  equation 

(20)  X*  -An*  +./A22 

for  L  is  opposite  in  stability  to  the  equation  (19)  for  K.  Thus  K(T  -  e) 
will  satisfy  the  Liapunov  equation 

(21)  An(T  -  e)K(T  -  e)  -  K(T  -  e)A22(T  -  e)  -  A12(T  -  e)  =  0  . 

A  A 

Its  solution  is  unique  because  A-j -j  and  A22  have  no  common  eigenvalues, 
cf.  e.g.  Bellman  [61.  An  explicit  solution  is  given  by  K(T  -  e)  = 

-M,2(T  -  e)Q22(T  -  e)  where  M^2  and  Q22  are  sub-blocks  of  the  modal 
matrix  M  for  A(T  -  e)  and  its  inverse  Q,  cf.  (11)  and  (14).  It  is 
preferable,  however,  to  obtain  K(T  -  e)  numerically  by  iteration  in  the 
equation 

(22)  Kj+1(T  -  e)  =  (An(T  -  e)K..(T  -  e)  -  A]2(T  -  e))A^(T  -  e) 

with  initial  iterate  K  (T  -  e)  =  0.  As  in  the  iteration  scheme  for 

o 

L(e),  the  convergence  will  be  rapid  for  |)A22(T  -  e)||  ||A-|^(T  -  e)[|  <_  y 


-  -/ 


small  in  the  spectral  norm.  When  A  is  time-invariant,  this  provides  the 
constant  matrix  K  appropriate  throughout  the  interval.  More  generally, 
however,  we  assume 


(H3)  The  solution  K(t)  of  the  linear  terminal  or  initial 
value  problem  (19),  (21)  stays  bounded  throughout 
0  <  t  <  T. 


We  note  that  (H3)  is  automatically  satisfied  if  K(T  -  e)  and  the  coeffi¬ 
cients  in  the  differential  equation  (19)  remain  bounded. 

With  these  three  hypotheses,  our  time-varying  LK  transformation  (5) 
has  now  become  completely  determined,  and  our  problem  (1)  is  reduced  to 
solving  the  time-scale  decoupled  system 


(23)  y}  =  A11(t)y1  +  B1(t)u  , 

(24)  y2  =  A22(t )y2  +  B2(t)u  , 

is  partitioned  after  its  first  n^  r 

Boundary  conditions  for  x  and  y  are  related  through  the  nonsingular 
matrix  T.  The  solutions  of  (23)  and  (24)  are  given  by 

t 


^1 

where  y  = 

and  B  =  TB  = 

-  2 

h 

ows. 


(25) 


yi(t)  =  Yi(t)ci  +  Yi(t)Yi:,(s)Bi(s)u(s)ds  , 


i  =  1  and  2,  for  constant  vectors  c^,  where  the  Y .  are  fundamental 
matrices  satisfying 


Y. 

i 


(0)  =  I 


ni 


Yi  =  A..(t)Y.  , 


0  <  t  <  T  . 
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Though  the  representation  (25)  is  useful,  direct  numerical  integration 
for  j/j  and  y2  Preferable  to  numerical  implementation  of  (25). 

Using  the  spectral  norm,  our  two-time-scale  assumption  implies  that 
jlA^U)!!  =  s(t)  while  1 1  A'^ft)  ( 1  £f_1(t),  so  j  |  A1 1  ( t)  j )  = 

!  i  A-j  -j  ( t )  I  j  I  I  ^22^22  i  1  1  s(t)f'](t))  jA22(t)  j  j  <  y(  |A22(t)|  |  .  Thus,  V2(t) 

/s 

is  rapidly  varying  compared  to  Y-j(t).  Indeed,  when  A22(t)  has  eigen¬ 
values  with  large  negative  real  parts,  say  of  order  0(~y),  Y2  decays  to 
zero  exponential ly  fast  and  it  becomes  negligible  outside  an  initial 

A 

0(cT)  boundary  layer.  Likewise,  when  the  eigenvalues  of  A^,(t)  are 
small,  like  0(k),  Yj(t)  is  nearly  constant  throughout  (0,T)  provided 
T  <<  1/k.  It  is  natural,  then,  to  think  of  y^  as  the  predominantly  slow 
solution  and  of  y^  as  the  predominantly  fast  solution,  realizing  that  the 
slow/fast  interpretation  could  be  corrupted  by  the  forcing  control 
B(t)u(t).  This  slow/fast  decomposition  would  carry  back  to  the  original 
system  (1)  as 


(26) 


x  ( t )  =  r>(t) 


ry1  (t)' 
y2(  t } 


Altogether,  then,  we've  transformed  our  original  problem  (1)  under 
hypotheses  (HI ) -( H3 )  into  the  integration  of  four  separate  problems  for 
L,  K,  y.j  and  y^,  with  L  and  K  being  constant  for  time-invariant  A 
matrices.  We'll  now  show  how  the  procedure  can  be  substantially 
simplified  through  approximations  when  we  impose  a  fast  mode  stability 
assumption.  Other  approximations  will  be  appropriate  under  different 
hypotheses. 
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3 .  Reduced  order  model inq  for  the  ini tial  value  problem— approximate 
analysi s. 

Let's  now  consider  the  initial  value  problem  for  (1),  assuming 
that  the  time  scale  parameter  (cf.  (3)  and  (HI))  is  small,  i.e. 

(27)  p  «  1  , 


and  that 


(H4)  the  eigenvalues  fj(t)  of  A^t)  all  have  large  negative 
real  parts  throughout  0  <  t  <  T. 


l<j<n2 

0<t<T 


Re  fj(t) 


also  holds. 


Because  (3)  implies  that  IjA^H  «  |  | A?2 1  [  ,  we  can  expect  the 

solution  of  the  linear  variational  equation  (20)  for  L  to  be  well- 

•  — 

approximated  through  the  nearby  system  L  =  -LA^.  Further,  the  large 

A 

magnitude  and  stability  of  the  eigenvalues  of  A^  suggest,  via  singular 
perturbations  theory,  that  the  initial  value  problem  for ^  will  have 
bounded  solutions  asymptotic  toJ£(t)  =  0  away  from  t  =  0,  while  the 
solution  of  the  corresponding  terminal  value  problem  will  become 
unbounded  for  t  <  T.  Therefore,  errors  made  in  the  numerical  integra¬ 
tion  of  the  Riccati  equation  for  L(t)  should  decay  exponential ly  to 
zero  in  forward  time  and  grow  exponentially  in  reverse  time.  To  keep 
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the  calculated  L(t)  bounded,  then,  under  hypothesis  (H4),  we  must  take 
e  =  0,  i.e.  we  define  L  and  K  through  initial  and  terminal  value 
problems,  respectively.  Indeed,  the  linear  system  for  K  will  be  well- 
approximated  through  the  nearby  system  K  =  -KA^  -  since  p  «  1, 
and  as  e  ->•  0  (and  IjA^Jj  -*-<»)  the  limiting  solution  will  satisfy 
K  =  0  for  t  <  T.  Thus,  the  familiar  quasi-steady  state  approximation, 
consistent  with  our  terminal  condition  K(T)  =  0,  holds  asymptotical ly. 
For  this  reason,  we  rewrite  the  system  (19)  for  K  as 

(29)  K(t)  =  K1 (t)  +  S(K( t ) ) 
with  the  nonhoinogeneous  term 

K-j  (t)  =  -A12(t)A^(t) 
and  the  linear  operator 

S(K)  =  (AnK  -  <)A"J  . 

We  shall  solve  the  system  by  successive  approximations,  starting  with 
the  trivial  iterate  KQ(t)  0.  Thus,  we  successively  define  the 
approximants 

(30)  K  K  +  j  S'(K,)  ,  j  >  2  , 

J  I  o  _  1  1 

for  K  where  Sf'(K^)  =  S(S~\k^))  for  each  i  >  1  and  S° ( )  =  K-j .  In 
practice,  only  a  few  iterates  will  be  needed  because  S(K)  has  a  small 
nonn  due  to  the  sizesofp  and  c,  i.e.  of  }  ( -| }  |  jjA^j]  and  JJA^JI. 
This  iteration  scheme  avoids  the  need  to  directly  integrate  the  terminal 


t 
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value  problem  for  K  and  to  store  its  solution  for  later  use  in  evaluating 
T  and  and  for  integrating  the  initial  value  problem  for  y^.  The 

successive  differentiations  of  involved  don't  pose  a  real  problem 
because  K(t)  is  asymptotically  negligible  for  t  <  T.  Indeed,  if  we  omit 
the  derivative  term  in  (29),  our  iterates  (30)  at  t  -  T  coincide  with 
those  of  (22)  used  to  obtain  K(T).  The  resulting  slow-mode  or  quasi¬ 
steady  state  approximation  K  (t)  to  K(t)  will  be  asymptotical ly  valid 
for  t  <  T.  The  approximation  K  (t)  -  K(t)  should  even  be  fairly  good 
near  t  =  T,  because  we  picked  K(T)  =  0. 

Returning,  then,  to  the  initial  value  problem  (24)  for  y^»  with 
y2 (0)  =[l(0)  I n x ( 0 ) ,  the  fact  that  Re  A(A22)  has  only  large  stable 
elements  suggests  that  y2  should  be  nearly  equal  to  its  slow-mode  (or 
quasi-steady  state)  approximation  y2  ;  0,  i.e. 

(31)  y2s(t)  -  -A^(t)B2(t)u(t) 

for  t  >  0.  Indeed,  singular  perturbations  theory  would  show  that  the 
"composite"  solution 

(32)  y2(t)  =  y2s(t)  +  y2f(t) 

will  provide  a  uniformly  valid  approximation  to  y2-  with  the  fast- 
varying  vector  y2^  satisfying 

(33)  y2f  =  A22(t)y2f  ,  y2f(0)  =  y2(0)  -  y2s(0) 

and  decaying  to  zero  exponentially  in  an  0(cT)  neighborhood  of  t  =  0. 

[If  a  good  approximation  to  y2  is  needed  near  t  =  0,  it  is  necessary 
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to  integrate  the  system  for  only  over  a  short  initial  interval,  but 
with  a  relatively  small  mesh  spacing. ]  Because  y^  satisfies 

•  /\  ^  ^ 

y2  =  A22y2  +  B2u  +  y2s  5  y2^  =  y2^ 

comparison  with  (24)  suggests  that  the  composite  vector  y^  will  be  a  good 
approximation  to  y2  (and  y^s  will  be  a  good  approximation  to  y^  away  from 
t  =  0)  provided  |  y2  |  is  small  on  0  £  t  <_  T  compared,  say,  to  the 
supretnum  of  { y2 ( 0) |  and  |y2s(t)j.  Thus,  we'll  assume 

(H5)  the  slow-mode  approximation  y^s  _to  y^  is  si owly-varyi ng 
throughout  0  <_  t  <_  T. 

We  recall  that  slowly-varying  functions  play  an  important  role  in 

asymptotic  analysis  (cf.  Feshchenko  et  ah  [13])  and  note  that  the 

assumption  is  reasonable  in  the  common  situation  that  y^s  is  itself 

small  when  |B2uj  is  small  compared  to  the  large  fjA^jj'.  Hypothesis  (H5) 

also  reflects  the  fact  that  rapid  variation  of  B^u  could  cause  the  state 

y2  to  be  fast  for  t  >  0,  even  though  the  free  response  would  be 

asymptotically  negligible  there.  We  note,  in  particular,  that  because  B^ 

[L  I  ]B,  y„  could  become  rapidly  varying  when  our  asymptotics  break- 
n  ^  cS 

down  because  L  isn't  small  or  the  forcing  Bu  is  rapid.  The  asymptotic 
decomposition  (32)  of  y2  into  slow  and  fast  parts  could  also  be  motivated 
by  using  Laplace's  method  (cf.  Olver  [25])  on  the  integral  representation 
(25)  for  y2- 

We  shall  integrate  the  full  n-j  dimensional  system  (23)  for  y^  using 

the  initial  vector  y,(0)  =  (I  +  K(0)L(0)  K(0))x(0).  (if  the  eigen- 

‘  1 

values  of  have  large  real  parts,  we  might  also  be  able  to  approximate 


* 
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y-j  by  a  quasi-steady  state  approximation  y.j  within  (0,T).  Then, 
however,  a  change  of  time  scale  s  =  At  for  an  appropriate  constant  A 
might  eliminate  this  stiffness.)  When  j|A^j|  isn't  large  and  the 
control  term  B^u  isn't  rapidly-varying, (23)  can  be  integrated  with  step- 
sizes  much  larcer  than  would  be  necessary  for  integrating  the  original 
system. 

For  t  >  0,  then,  the  solution  of  our  original  problem  will  be  well- 
approximated  by  the  slow-mode  approximation 

(t) " 

(34>  J(t)  ■  r'’l(t)|y2s<t)  ■ 

(If  desired,  we'd  have  to  correct  this  near  t  =  0  by  taking  the  fast¬ 
mode  correction  y^t)  to  y^s  into  account.)  For  t  >  0,  we've  achieved 
a  substantial  order  reduction  because  we  need  only  integrate  an  initial 
value  problem  for  l(t)  and  another  for  y-j(t).  This  is  because  is 
obtained  explicitly  from  the  algebraic  equation  (31)  and  K  is  obtained 
from  a  fast-converging  iteration  scheme  (30),  under  our  fast-mode 
stability  assumption. 

All  these  arguments  can  be  made  completely  rigorous  by  explicitly 
using  the  small  parameters  u  and  c  to  rescale  our  differential  equations 
and  carrying  out  a  careful  asymptotic  analysis  as  c  and  p  simultaneously 
tend  toward  zero.  For  only  moderately  small  parameters,  a  full  integra¬ 
tion  of  the  linear  systems  for  K  and  y->  might  be  needed.  The  connection 
between  our  approximations  and  numerical  methods  for  systems  of  stiff 
differential  equations  is  closest  to  the  smooth  approximate  particular 
solution  techniqueof  Dahlquist  [9]  and  Oden  [24]. 


.. 


To  summarize,  we  list  the  somewhat  oversimplified  steps  appropriate 


for  obtaining  a  reduced  order  model  of  our  fast-mode  stable,  two-time 
scale  system  on  t  >  0.  They  are: 


(1)  Use  the  eigenvalues  of  A(0)  to  determine  the  number  n^  of 
slow  modes. 

(2)  Obtain  L(0)  by  iterating  in  the  equation 

Li  +  1  =  (A22(0)  +  LiA12(0))'1(LiAn(0)  +  A21(0)), 

i  >  0,  Lq  -  0. 

(3)  Integrate  the  initial  value  problem  for 

L  =  A^L  -  LA^-j  +  LA^L  -  A^,  making  sure  that  the  trans¬ 
formed  system  remains  fast-mode  stable  and  two-time  scale 
with  n^  constant  throughout  0  <_  t  £_  T  and  that  the  slow¬ 
mode  approximation  y^s  remains  slowly-varying  there. 

(4)  Obtain  K(t)  on  [0,T]  through  the  iteration 
Kj+1 ( t) =  ('A,2(t)  +  A1 1 ( t)Kj ( t)  -  Kj(t))A'2(t), 

j  >  0,  K  ( t )  =  0.  (Alternatively,  obtain  the  slow-mode 
approximation  K  (t)  by  omitting  the  derivative  term.) 

(5)  Integrate  the  initial  value  problem  for  y^(t)  and  obtain  the 
reduced-order  solution 


x  ( t ) 


'  -  L  ( t ) 


lyl(t)  + 


-K(t) 
h  L(t)K(t)J 


!yzs ( t)  for  t  > 
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4.  Related  Problems. 

a.  Nearly  constant  slow  modes  are  found  for  two-time-scale  systems 

when  the  eigenvalues  of  are  not  large.  Then,  the  small  size  of  the 

eigenvalues  of  suggests  that  Y^(t)  is  nearly  constant  on  a  fixed 

finite  interval.  Though  the  dynamics  for  y^  are  not  simplified,  we 

obtain  order  reduction  in  the  sense  that  y-j  will  simply  track  its 

t 

forcing,  i.e.  y^(t)  =  y-j(O)  +  j  B-j  (s)u(s)ds . 

0 

b.  It  may  sometimes  be  simpler  to  simply  block-triangularize  our 
system  matrix  through  the  matrix  T  in  the  system  (9)  for 


z  =  ,  the  fast  modes  are  decoupled  via  L(t),  so  a  slow-mode  approxi- 

mation  for  z^  could  be  used  in  the  forcing  for  z^  on  t  >  0.  We  block- 
diagonalized  our  system,  since  the  linear  problem  for  K(t)  seems  simple 
after  the  quadratic  problem  for  L(t). 

c.  When  the  eigenvalues  of  A^  have  both  large  positive  and  large 
negative  real  parts,  the  initial  (terminal)  value  problem  for  L(t)(K(t)) 
will  no  longer  be  well -posed.  Only  certain  two-point  problems  for  x  can 
be  expected  to  have  bounded  solutions  (cf.  O'Malley  [26]  and  O'Malley  and 
Anderson  [28]).  Effective  use  of  time-scale  separation  should,  nonethe¬ 
less,  be  computationally  significant  in  obtaining  approximate  solutions 
to  appropriate  two-point  problems. 
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d.  The  n-|  "slow"  solutions  of  the  unforced  problem  are  spanned  by 

‘nr 

the  columns  of  the  matrix  Y,(t).  If  we  therefore  integrate  the 

[-L(t)j  1 

initially  slow  modes  of  our  system  (1)  forward  in  time  to  obtain  the 

X(t)'i 

j,  we'll  necessarily  have  L(t)  =  -X0(t)X~  (t).  Thus 
(X2(t)j  2  1 


n  x  n-j  matrix 


existence  of  L  is  guaranteed  as  long  as  X-j  ( t )  remains  nonsingular.  For 
problems  where  L  becomes  unbounded  on  0  <  t  <  T,  there  still  remains  the 
possibility  of  reinitializing  our  problem  to  keep  the  appropriate  n-j  *  n^ 
matrix  nonsingular.  This  corresponds  to  the  reorthonormalizations  used 
by  Scott  and  Watts  [32] . 


5 .  Numerical  Examples. 

In  practice,  the  need  for  reduced-order  modelling  requires  us  to  use 
our  schemes  on  problems  where  the  time-scale  separation  parameter  u  and 
the  fast-mode  stability  parameter  e  are  not  asymptotically  small.  Among 
many  other  considerations,  we  must  then  be  particularly  concerned  with 
the  difference  between  eigenvalue  stability  and  actual  stability  and  with 
the  occurrence  of  eigenvalues  with  large  imaginary  parts  that  can  allow 
slow  modes  (so  classified  by  eigenvalue  magnitudes)  to  decay  faster  than 
some  fast-modes.  The  latter  concern  might  be  illustrated  through  a  third 
order  system  with  the  slow  eigenvalues  s  =  -1  and  the  fast  eigenvalues 
f.j  ^  =  -0.1  +  ilO.  Then,  c  =  lj?-,  so  for  T  sufficiently  large,  all  modes 
will  be  negligible  away  from  t  =  0.  Otherwise,  the  fast  modes  cannot  be 
ignored.  A  check  on  the  relative  stability  of  the  slow  and  fast 


subsystems  can  be  made  through  the  ratio  o(t)  =  -max j Re(s . ( t) ) | 

i 

/max(Re  f - ( t) )  and  its  maximum  over  0  <  t  <  T. 

J  J  '  , 

If  A(t)  has  the  time-varying  spectral  decomposition  A  =  MJM  ,  the 

change  of  variables  w  =  M_1x  converts  the  problem  (1)  into  w  =  (J  -  M_1M)w 

M  ]8u.  Thus  eigenvalue  rotation,  measured  by  the  size  of  M_1M,  can 

substantially  alter  the  stability  suggested  by  the  eigenvalues  of  A  and  J. 

Rapid  variation  of  the  slow-eigenspace  of  A  could,  in  particular,  make  L 

and  y2$  rapidly-varying,  and  jeopardize  the  appropriateness  of  our 

approximations . 

We  shall  consider  two  time- varying  third  order  problems  with  one 
slow  mode.  Specifically,  let  the  state  matrices  A.  =  MJ.M”1  for  i  -  1 
and  2  have 


A 

+  h(t) 

0 

0  1 

0 

1 

0  1 

0 

0 

(1  +  h(t))'1] 

with  and  being  real  canonical  forms  with  snectra  X(J-|)  = 

-(1  +  h 1  ( t) ) {1 , 10 ,1 2)  and  A ( J 2 )  =  -(1  +  h'(t))fl,3  +  10i}.  Thus  for  both 

examples,  p  «  0.1  while  ( 1  -j , r 2 )  =  (0-022,0.074)  for  T  =  4.  Therefore,  the 

fast-mode  stability  and  the  relative  stability  of  the  fast  modes  is  stronger 

'll  o 

for  the  first  example.  We'll  also  take  x(0)  =  3  =  1  ,  u(t)  =  sin  irt, 

, 

and  h(t)  =  h'(t)  =  <1-  sin-;t/4.  ('J 

O 

For  both  examples,  the  appropriate  initial  condition  for  the  Riccati 

ft,  (on  fn 

differential  equation  (10)  is  the  two  vector  L(0)  =  ,ni  =  I  1 .  Since 

i2! 

the  quadratic  equation  (12)  provides  a  steady-state  for  the  corresponding 
differential  equation  (10)  at  t  =  0,  we  might  attempt  to  find  L(0)  as  an 
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equilibrium  solution.  Figure  5  reoresents  the  phase  plane  for 

example  one.  As  shown,  all  points  above  the  line  >’ ^  =  3^  -  1  converge 
to  L(0),  but  points  below  this  line  diverge  to  infinity.  For  example 
two,  a  slow  oscillatory  convergence  is  illustrated  in  Figure  6.  Thus, 
this  natural  way  to  seek  L ( 0)  is  only  locally  convergent,  in  contrast  to 
the  safer,  globally  convergent  iterative  method  we  described  previously 
via  (15). 

Once  L ( 0 )  is  obtained,  the  time-varying  Riccati  equation  (10)  can 

ft,(t)' 

be  integrated  from  t  =  0  to  4.  The  solution  L(t)  =  L  for  example 

V'2 

one  is  illustrated  in  Figure  7.  We  have  also  plotted  the  smooth  solution 
L ( t)  with  L ( 0 )  =  L(0)  of  the  algebraic  Riccati  equation  obtained  when  we 
set  the  derivative  term  in  (10)  to  zero.  L(t)  is  a  good  approximation  to 
L(t).  This  should  not  be  unexpected  since  E(t)  =  L(t)  -  L(t)  satisifes 
E  =  (A^2  +  LA^)E  -  E  ( A  -j  -j  -  A-^L)  +  EA-^E  -  L  on  0  <  t  <  4  with  E(0)  =  0. 
Presuming  A  22  +  LA^  maintains  large,  strongly  stable  eigenvalues 
compared  to  A^  -  A^L  and  presuming  L  isn't  large,  singular  perturbations 
would  suggest  a  small  error  E(t)  throughout  0  <_  t  <_  4.  Thus,  we  could 
often  expect  to  use  L,  the  solution  of  an  algebraic  system,  to  approximate 

A 

L(t).  Figure  7  also  includes  trajectories  L(t)  for  the  Riccati  system 
with  perturbed  initial  conditions.  They,  too,  converge  to  L(t)  for  t  >  0 
provided  the  initial  perturbations  aren't  too  large.  With  the  weaker 
fast-mode  stability  of  Example  2,  we  found  that  the  initial  behavior  of  L 
trajectories  was  oscillatory  and  convergence  to  L  was  delayed  (cf. 

Figure  8). 


- 


Figure  5: 


The  phase  plane 
for  the  Riccati 
solution  component? 
for  Example  1. 


Figure  6: 

The  phase  plane 
for  the  Riccati 
solution  component 
for  Example  2. 


Figure  7: 


Solutions  l(4 
and  L(t)  of 
the  Riccati 
differential 
equation  and 
the  algebraic 
Riccati 
equation  for 
Example  1 . 


Figure  8: 

Solutions  L( 
and  L ( t )  of 
the  Riccati 
differential 
equation  anc 
the  algebrai 
Riccati 
equation  for 
Example  2. 


The  linear  differential  system  (19)  for  K(t)  =  ( k ^ ( t )  ( t ) )  was 

integrated  backward  in  time  from  t  =  4  to  t  =  0.  It  could  also  be 
solved  readily  via  the  iteration  approach.  The  relative  behavior  of 
K(t),  of  the  corresponding  linear  algebraic  problem  for  K(t),  and  of  the 
problem  for  K(t)  with  perturbed  terminal  values  is  analogous  to  that  for 

^  A  A 

L(t),  L(t),  and  t-(t),  except  that  the  convergence  of  K  to  K  also  holds 
for  large  perturbations  of  end  vectors. 

By  changing  h'  to  the  more  oscillatory  g  sin  nt,  the  eigenvalues  of 
A  are  changed,  but  there  is  little  change  in  the  solution  L(t)  of  the 
Riccati  system.  The  corresponding  change  of  h(t)  to  ^  sin  nt,  however, 
produces  a  more  rapid  oscillation  in  the  eigenspace  of  A(t)  and  there 
results  more  rapid  change  in  the  decoupling  vector  L(t).  Nonetheless, 
as  Figures  9  and  10  illustrate,  L(t)  still  remains  a  good  approximation 
to  L(t) . 

By  superposition,  the  solution  x(t)  of  our  forced  initial  value 
problem  (1)  can  be  considered  to  be  the  sum  of  the  separate  responses  of 
the  unforced  system  with  initial  vector  x(0)  and  to  the  input  u(t)  with 
zero  initial  state.  The  unforced  response  is  illustrated  in  Figures  11 
and  12.  The  exact  solution  x(t)  is  well  approximated  by  the  first  order 
approximation  x(t)  outside  an  initial  transient  region  of  approximate 
thickness  20c.  This  corresponds  to  0.5  time  units  for  Example  1  and  1.7 
units  for  Example  2.  For  slowly-varying  control  inputs  u(t),  the 
agreement  between  x(t)  and  x(t)  is  very  good  for  both  examples.  For  the 
rapidly  varying  control  input  u(t)  =  sin  ut,  however,  the  resulting 
approximations  x(t)  for  Example  1  are  better  than  for  Example  2,  which 


has  weaker  fast-mode  stability  (cf.  Figures  13  and  14). 


Figure  9: 


Riccati 
solutions  foi 
Example  1  am 
a  more 
osci 1 latory 
h'(t). 


Figure  10: 

Riccati  solu¬ 
tions  for 
Example  2  and 
new  h1 (t) . 


Figure  11 


Response  of 
Example  1 
with  u(t)  = 


Figure  12: 

Response  of 
Example  2  wi 
u(t)  =  0. 


time,  t 


Figure  13: 


ResDonse  of 
Example  1  to 
control 

2 

u ( t )  =  sin  irt 
and  initial 
state  x(0)  =  0 


time,  t 


Figure  14: 

Response  of 
Example  2  to 
control 

2 

u(t)  =  sin  n 
and  initial 
state  x(0)  - 


3.0 


2.0 


time,  t 
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